function [W_Gi,params] = temporal_sphericalwave(params,Z)
% params: parameters
% V: matrix comes from SVD
% Z: matrix for saving j, n and m

lambda = params.lambda;
W_V = params.W_V;
Cn2 = params.Cn2;
L = params.L;
D = params.D;
t_miu = params.t_miu;
f_y = params.f_y;
k = 2*pi/lambda;
layerNum = params.layerNum;
L0 = params.L0; % Outer scale
l0 = params.l0; % Inner scale

H = (0:layerNum)*L/layerNum;
H_bot = H(1:end-1);
H_top = H(2:end);
% dimensions: [m], [m^(1/3)], [m/s]
MLM = zeros(layerNum,4);
MLM(:,1) = (H_top-H_bot)'; % 分块区间的区间大小 [m]
for i = 1:layerNum
    H_mid = H_bot(i) + (H_top(i)-H_bot(i))/2;
    if (strcmp(class(Cn2),'double'))
        MLM(i,2) = Cn2;
        MLM(i,3) = W_V;
    else
        intFun1 = @(h) Cn2(h*10^(-3));
        intFun2 = @(h) W_V(h*10^(-3));
        MLM(i,2) = 1/MLM(i,1)*integral(intFun1,H_bot(i),H_top(i)); %Cn2_q
        MLM(i,3) = 1/MLM(i,1)*integral(intFun2,H_bot(i),H_top(i)); %V_q
    end
    MLM(i,4) = H_mid; %z_q
end
params.MLM = MLM;

% multilayer model for r0 parameter
r0m = zeros(1,layerNum);
for i = 1:layerNum
    r0m(i) = r0_generate(Cn2,lambda,H_bot(i),H_top(i));
end
params.r0m = r0m;
r0 = r0m_merge(r0m);
params.r0 = r0;

% computing M_i
W_Gi = 0;
intFun3 = @(h) (h/L).^(8/3);
for j = 1:layerNum
    M_i = K_L_to_Mi(Z, t_miu/MLM(j,3)*(MLM(j,4)/L), f_y, D);
    [t_miu2, f_y2] = meshgrid(t_miu/MLM(j,3)*(MLM(j,4)/L), f_y);
    fun_cal = (abs(M_i).^(2)).*(((t_miu2).^2+(f_y2).^2+L0^(-2)).^(-11/6))...
            .*exp(-(2*pi*l0/5.92)^2*((t_miu2).^2+(f_y2).^2));
    cal_integral = sum(fun_cal,1)*(1/D);
    F_Gi = 0.033*(2*pi)^(-2/3)*(2*pi/lambda)^2*cal_integral;
    W_Gi = (MLM(j,2)/MLM(j,3))*F_Gi...
        *integral(intFun3,H_bot(j),H_top(j)) + W_Gi;
    fprintf('multilayer calc %3g / %3g \n', j, layerNum);
end

% end of the function
end